{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Axisymmetric Newtonian incompressible Navier-Stokes element formulation and automatic differentiation\n",
    "**Rubén Zorrilla**\n",
    "UPC BarcelonaTech, CIMNE\n",
    "\n",
    "## Governing equations\n",
    "Starting from the expression of the incompressible Navier-Stokes equations for a Newtonian fluid in polar coordinates and applying the simplifications corresponding to an axisymmetric flow one obtains the linear momentum conservation equations\n",
    "\n",
    "\\begin{align*}\n",
    "\\rho\\frac{\\partial u_{x}}{\\partial t} + \\rho u_{y}\\frac{\\partial u_{x}}{\\partial y} + \\rho u_{x}\\frac{\\partial u_{x}}{\\partial x} = & - \\frac{\\partial p}{\\partial x} + \\mu\\left(\\frac{1}{y}\\frac{\\partial}{\\partial y}\\left(y\\frac{\\partial u_{x}}{\\partial y}\\right) + \\frac{\\partial^{2}u_{x}}{\\partial x^{2}}\\right) + \\rho g_{x} \\; , \\\\\n",
    "\\rho\\frac{\\partial u_{y}}{\\partial t} + \\rho u_{y}\\frac{\\partial u_{y}}{\\partial y} + \\rho u_{x}\\frac{\\partial u_{y}}{\\partial x} = & - \\frac{\\partial p}{\\partial y} + \\mu\\left(-\\frac{u_{y}}{y^{2}} + \\frac{1}{y}\\frac{\\partial}{\\partial y}\\left(y\\frac{\\partial u_{y}}{\\partial y}\\right) + \\frac{\\partial^{2}u_{y}}{\\partial x^{2}}\\right) + \\rho g_{y}\n",
    "\\end{align*}\n",
    "and the mass conservation equation\n",
    "\\begin{equation*}\n",
    "\\frac{1}{y}\\frac{\\partial}{\\partial y}\\left(y u_{y}\\right) + \\frac{\\partial u_{x}}{\\partial x} = 0\\; .\n",
    "\\end{equation*}\n",
    "\n",
    "Note that this implementation assumes x-coordinate to be aligned with the revolution axis, meaning that y-coordinate represents the radial velocity.\n",
    "Also note that current derivation neglects the viscous isochoric stress coming from the incompressibility error.\n",
    "\n",
    "## Variational form\n",
    "The variational form (i.e. the functional to be automatically differenctiated) is obtained by multiplying previous equations by the corresponding test functions $\\mathbf{w}$ and $q$. After integration by parts, this yields\n",
    "\n",
    "\\begin{align*}\n",
    "\\left(w_{x},\\rho g_{x}\\right)_{\\Omega} - & \\left(w_{x},\\rho\\frac{\\partial u_{x}}{\\partial t}\\right)_{\\Omega} - \\left(w_{x},\\rho a_{y}\\frac{\\partial u_{x}}{\\partial y}\\right)_{\\Omega} - \\left(w_{x},\\rho a_{x}\\frac{\\partial u_{x}}{\\partial x}\\right)_{\\Omega} + \\left(\\frac{\\partial w_{x}}{\\partial x}, p\\right)_{\\Omega} \\\\\n",
    "- & \\left(\\frac{\\partial w_{x}}{\\partial y}, \\mu\\frac{\\partial u_{x}}{\\partial y}\\right)_{\\Omega} - \\left(\\frac{\\partial w_{x}}{\\partial x}, \\mu\\frac{\\partial u_{x}}{\\partial x}\\right)_{\\Omega} + \\langle w_{x}, \\mu\\left(\\frac{\\partial u_{x}}{\\partial y}+\\frac{\\partial u_{x}}{\\partial x}\\right)n_{x} - p n_{x} \\rangle_{\\Gamma} \\\\\n",
    "+ & \\left(w_{y},\\rho g_{y}\\right)_{\\Omega} - \\left(w_{y},\\rho\\frac{\\partial u_{y}}{\\partial t}\\right)_{\\Omega} - \\left(w_{y},\\rho a_{y}\\frac{\\partial u_{y}}{\\partial y}\\right)_{\\Omega} - \\left(w_{y},\\rho a_{x}\\frac{\\partial u_{y}}{\\partial x}\\right)_{\\Omega}\n",
    "+ \\left(\\frac{w_{y}}{y} + \\frac{\\partial w_{y}}{\\partial y}, p\\right)_{\\Omega} \\\\\n",
    "- & \\left(\\frac{\\partial w_{y}}{\\partial y}, \\mu\\frac{\\partial u_{y}}{\\partial y}\\right)_{\\Omega} - \\left(\\frac{\\partial w_{y}}{\\partial x}, \\mu\\frac{\\partial u_{y}}{\\partial x}\\right)_{\\Omega} + \\langle w_{y}, \\mu\\left(\\frac{\\partial u_{y}}{\\partial y}+\\frac{\\partial u_{y}}{\\partial x}\\right)n_{y} - p n_{y} \\rangle_{\\Gamma} \\\\\n",
    "- & \\left(q,\\frac{u_{y}}{y}\\right)_{\\Omega} - \\left(q,\\frac{\\partial u_{y}}{\\partial y}\\right)_{\\Omega} - \\left(q, \\frac{\\partial u_{x}}{\\partial x}\\right)_{\\Omega} = 0\\; ,\n",
    "\\end{align*}\n",
    "being $\\mathbf{a}$ the convective velocity, which computed from the linearised previous iteration velocity $\\mathbf{u}^{k-1}$ and the mesh velocity $\\mathbf{u}_{M}$ as\n",
    "\\begin{equation*}\n",
    "\\mathbf{a} = \\mathbf{u}^{k-1} - \\mathbf{u}_{M}\\; .\n",
    "\\end{equation*}\n",
    "\n",
    "Note that the cylindrical component of the divergence operator applied to the viscous stress term (see the strong form above) no longer appears in the variational form as we get rid of it by integrating by parts.\n",
    "However, it appears in the test function divergence of the pressure term.\n",
    "\n",
    "## Variational MultiScales stabilization\n",
    "The scales separation reads\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathbf{u} & = \\mathbf{u_{h}}+\\mathbf{u_{s}}\\\\\n",
    "p & = p_{h}+p_{s}\n",
    "\\end{align*}\n",
    "\n",
    "being $\\mathbf{u}_{h}$ and $p_{h}$ the finite element scales and $\\mathbf{u}_{s}$ and $p_{s}$ the subgrid scales, which are modelled from the finite element residuals as\n",
    "\n",
    "\\begin{align*}\n",
    "u_{s,x} & = \\tau_{1} R^{M}_{x}(u_{h,x},p_{h})\\; , \\\\\n",
    "u_{s,y} & = \\tau_{1} R^{M}_{y}(u_{h,y},p_{h})\\; , \\\\\n",
    "p_{s} & = \\tau_{2} R^{C}(\\mathbf{u}_{h},p_{h})\\; ,\n",
    "\\end{align*}\n",
    "\n",
    "being $\\tau_{1}$ and $\\tau_{2}$ the stabilization parameters which in this case equal\n",
    "\\begin{equation*}\n",
    "    \\tau_{1} = \\left(\\frac{\\rho\\tau_{d}}{\\Delta t} + \\frac{c_{1}\\mu}{h^{2}} + \\frac{c_{2}\\rho\\lVert\\mathbf{a}\\rVert}{h}\\right)^{-1}\n",
    "\\end{equation*}\n",
    "and\n",
    "\\begin{equation*}\n",
    "    \\tau_{2} = \\mu+\\frac{c_{2}\\rho h\\lVert\\mathbf{a}\\rVert}{c_{1}} \\; ,\n",
    "\\end{equation*}\n",
    "being $c_{1}$ and $c_{2}$ the user-definable stabilization constants, $\\tau_{d}$ a dynamic coefficient ranging between 0 and 1 and $h$ the characteristic element size.\n",
    "\n",
    "The finite element residuals $\\mathbf{R}^{M}(\\mathbf{u}_{h},p_{h}$ and $R^{C}(\\mathbf{u}_{h},p_{h})$ are defined as\n",
    "\\begin{align*}\n",
    "R^{M}_{x}(u_{h,x},p_{h}) = & \\rho g_{x} - \\rho\\frac{\\partial u_{h,x}}{\\partial t} - \\rho a_{h,y}\\frac{\\partial u_{h,x}}{\\partial y} - \\rho a_{h,x}\\frac{\\partial u_{h,x}}{\\partial x} - \\frac{\\partial p_{h}}{\\partial x} + \\mu\\left(\\frac{1}{y}\\frac{\\partial u_{h,x}}{\\partial y} + \\frac{\\partial^{2} u_{h,x}}{\\partial y^{2}} + \\frac{\\partial^{2} u_{h,x}}{\\partial x^{2}}\\right)\\; , \\\\\n",
    "R^{M}_{y}(u_{h,y},p_{h}) = & \\rho g_{y} - \\rho\\frac{\\partial u_{h,y}}{\\partial t} - \\rho a_{h,y}\\frac{\\partial u_{h,y}}{\\partial y} - \\rho a_{h,x}\\frac{\\partial u_{h,y}}{\\partial x} - \\frac{\\partial p_{h}}{\\partial y} + \\mu\\left(- \\frac{u_{h,y}}{y^{2}} + \\frac{1}{y}\\frac{\\partial u_{h,y}}{\\partial y} + \\frac{\\partial^{2} u_{h,y}}{\\partial y^{2}} + \\frac{\\partial^{2} u_{h,y}}{\\partial x^{2}}\\right)\\; , \\\\\n",
    "R^{C}(\\mathbf{u}_{h},p_{h}) = & -\\frac{1}{y}u_{h,y} - \\frac{\\partial u_{h,y}}{\\partial y} - \\frac{\\partial u_{h,x}}{\\partial x} \\; .\n",
    "\\end{align*}\n",
    "\n",
    "Once the subscales are defined, one can derive the stabilization terms to be added to the functional above. Hence, after doing the required integration by parts to remove the subscale derivatives and dropping the high order terms (the use of linear elements is assumed), one obtains the stabilization functional\n",
    "\\begin{align*}\n",
    "\\left(\\frac{\\partial w_{x}}{\\partial y}\\rho a_{y} + w_{x}\\rho\\frac{\\partial a_{y}}{\\partial y}, u_{s,x}\\right)_{\\Omega} & + \\left(\\frac{\\partial w_{x}}{\\partial x}\\rho a_{x} + w_{x}\\rho\\frac{\\partial a_{x}}{\\partial x}, u_{s,x}\\right)_{\\Omega} + \\left(\\frac{\\partial w_{x}}{\\partial x},p_{s}\\right)_{\\Omega} \\\\\n",
    "& + \\left(\\frac{\\partial w_{y}}{\\partial y}\\rho a_{y} + w_{y}\\rho\\frac{\\partial a_{y}}{\\partial y}, u_{s,y}\\right)_{\\Omega} + \\left(\\frac{\\partial w_{y}}{\\partial x}\\rho a_{x} + w_{y}\\rho\\frac{\\partial a_{x}}{\\partial x}, u_{s,y}\\right)_{\\Omega} + \\left(\\frac{w_{y}}{y} + \\frac{\\partial w_{y}}{\\partial y},p_{s}\\right)_{\\Omega} \\\\\n",
    "& - \\left(q,\\frac{1}{y}u_{s,y}\\right)_{\\Omega} + \\left(\\frac{\\partial q}{\\partial y},u_{s,y}\\right)_{\\Omega} + \\left(\\frac{\\partial q}{\\partial x},u_{s,x}\\right)_{\\Omega}\n",
    "\\end{align*}\n",
    "in which the quasi-static subscales assumption has been considered.\n",
    "\n",
    "## IMPLEMENTATION\n",
    "The remaining implementation is similar to that of the weakly compressible Navier-Stokes element but with the added simplification of assuming a Newtonian constitutive model within the element. \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
